require(mapdata)
require(reshape)
rm(list=ls())
setwd('/media/Iomega_HDD/boulot/Article/kd/last_23_01_2012/data')  		##Choose your path

data<-'Global_Data_2011_12_14v3.txt'   ##name of the dataset file

###################################################################################
################################  get the dataset  ################################
###################################################################################
t<-read.table(data, header=T)

####################  Table has the following column names  #######################
#  "station_id"    "cruise_id"     "experiment"    "investigator"  "cruise"        "date"          "year"
#  "month"         "day"           "start_time"    "end_time"      "min_latitude"  "max_latitude"  "min_longitude"
#  "max_longitude" "lat"           "lon"           "etopo2"        "chla_f"        "chla_h"        "chla_hf"      
#  "chla_fh"       "chla_avg"      "gsm_chl"       "gsm_cdm"       "gsm_bbp"       "Rrs320"        "Rrs340"       
#  "Rrs380"        "Rrs412"        "Rrs443"        "Rrs455"        "Rrs465"        "Rrs490"        "Rrs510"       
#  "Rrs520"        "Rrs530"        "Rrs550"        "Rrs555"        "Rrs560"        "Rrs565"        "Rrs570"       
#  "Rrs590"        "Rrs620"        "Rrs625"        "Rrs665"        "Rrs670"        "Rrs683"        "Rrs700"       
#  "rrs_qc"        "ad320"         "ad340"         "ad380"         "ad412"         "ad443"         "ad455"        
#  "ad465"         "ad490"         "ad510"         "ad520"         "ad530"         "ad550"         "ad555"        
#  "ad560"         "ad565"         "ad570"         "ad590"         "ad620"         "ad625"         "ad665"        
#  "ad670"         "ad683"         "ad700"         "ap320"         "ap340"         "ap380"         "ap412"        
#  "ap443"         "ap455"         "ap465"         "ap490"         "ap510"         "ap520"         "ap530"        
#  "ap550"         "ap555"         "ap560"         "ap565"         "ap570"         "ap590"         "ap620"        
#  "ap625"         "ap665"         "ap670"         "ap683"         "ap700"         "ap_qc"         "ag320"        
#  "ag340"         "ag380"         "ag412"         "ag443"         "ag455"         "ag465"         "ag490"        
#  "ag510"         "ag520"         "ag530"         "ag550"         "ag555"         "ag560"         "ag565"        
#  "ag570"         "ag590"         "ag620"         "ag625"         "ag664"         "ag670"         "ag683"        
#  "ag700"         "ag_qc"         "bb405"         "bb412"         "bb443"         "bb448"         "bb455"        
#  "bb465"         "bb470"         "bb490"         "bb510"         "bb520"         "bb530"         "bb550"        
#  "bb555"         "bb560"         "bb565"         "bb570"         "bb590"         "bb620"         "bb625"        
#  "bb650"         "bb660"         "bb665"         "bb670"         "bb676"         "bb683"         "bb715"        
#  "bb_qc"         "Kd320"         "Kd340"         "Kd380"         "Kd412"         "Kd443"         "Kd455"        
#  "Kd465"         "Kd490"         "Kd510"         "Kd520"         "Kd530"         "Kd550"         "Kd555"        
#  "Kd560"         "Kd565"         "Kd570"         "Kd590"         "Kd620"         "Kd625"         "Kd665"        
#  "Kd670"         "Kd683"         "Kd700"         "kd_qc"         "cdm_slope"     "bbp5_expo"     "bbp7_expo"

### First, we clean the data to only keep the line where there is at least one value of Kd and GSM value

wavelength <- c(320,340,380,412,443,490,510,555)
kd_names <- paste("Kd", wavelength,sep="")
index_names <- which(colnames(t)%in%kd_names)
t2 <- t[,index_names]
t2[t2 < -5]<-NA

na_func <- function(t2){
  sum(is.na(t2))  
}
nbna<-apply(t2,1,na_func)

## we just keep the lines where there is at least a Kd value (for any wavelength) and chl are positive
t3 <- t[-which(nbna==length(wavelength) & t$gsm_chl<0 & t$chla_avg<0  &
  t$kd_qc!=0),]
save(t3, file="t3.Rdata")

t4<-cbind(t3$cruise,t3$month,t3$day,t3$year,t3$lat,t3$lon,t3$chla_avg,t3$gsm_chl,t3$gsm_cdm,
          t3$gsm_bbp,t3$Kd320,t3$Kd340,t3$Kd380,t3$Kd412,t3$Kd443,t3$Kd490,t3$Kd510,t3$Kd555)
t4[t4==-999] <- NA 
colnames(t4) <- c("cruise", "month", "day", "year", "lat", "lon", "chla_avg", "gsm_chl", "gsm_cdm", "gsm_bbp",
                  "Kd320", "Kd340", "Kd380", "Kd412", "Kd443", "Kd490", "Kd510", "Kd555")
save(t4, file="t4.Rdata")

t5 <- melt(as.data.frame(t4), measure.vars=colnames(t4)[11:18])
t5 <- t5[-which(is.na(t5$value)==T | is.na(t5$gsm_chl)==T),]

t5$lambda <- 1
t5$convert_lambda <- 1
t5$lambda[which(t5$variable=="Kd320")] <- 320; t5$convert_lambda[which(t5$variable=="Kd320")] <- 1
t5$lambda[which(t5$variable=="Kd340")] <- 340; t5$convert_lambda[which(t5$variable=="Kd340")] <- 2
t5$lambda[which(t5$variable=="Kd380")] <- 380; t5$convert_lambda[which(t5$variable=="Kd380")] <- 3
t5$lambda[which(t5$variable=="Kd412")] <- 412; t5$convert_lambda[which(t5$variable=="Kd412")] <- 4
t5$lambda[which(t5$variable=="Kd443")] <- 443; t5$convert_lambda[which(t5$variable=="Kd443")] <- 5
t5$lambda[which(t5$variable=="Kd490")] <- 490; t5$convert_lambda[which(t5$variable=="Kd490")] <- 6
t5$lambda[which(t5$variable=="Kd510")] <- 510; t5$convert_lambda[which(t5$variable=="Kd510")] <- 7
t5$lambda[which(t5$variable=="Kd555")] <- 555; t5$convert_lambda[which(t5$variable=="Kd555")] <- 8

save(t5, file="dataset.Rdata")

data_admb <- cbind(signif(t5$value,7), signif(t5$gsm_chl,7), signif(t5$gsm_cdm,7), signif(t5$gsm_bbp,7),
              t5$lambda, t5$convert_lambda)

##################################################################
###############Calculate the Kw ##################################
##################################################################
kw_MM3<-read.table('asw_bbsw_300_1100nm.txt',header=T,sep="\t")   ### from Morel & Maritorena (2001)

Kw <- kw_MM3$asw[which(kw_MM3$lambda%in%wavelength)] + kw_MM3$bbsw[which(kw_MM3$lambda%in%wavelength)]

### write the file for ADMB
filename <- "model.dat"

cat("###  init_number Nobs_tot                 //Total number of observations\n",
    dim(t5)[1],
    "\n###  init_vector Kw    //Value of Water absorption for each wavelength\n",
    file=filename, sep=""
)
write.table(paste(Kw, sep=" "), file=filename,append=T,eol=" ", row.names=FALSE, col.names=FALSE, quote=F,sep=' ')
cat("\n###  init_matrix Data(1,Nobs_tot,1,6)     //The data extracted using the R script\n",
    file=filename, sep="", append=T
)

write.table(data_admb,file=filename,append=T,eol="\n", row.names=FALSE, col.names=FALSE, quote=F,sep='\t')

file.copy(filename, paste("../script_admb/",filename, sep=""), overwrite=T)
#file.copy(fromfilename, to=paste("~/Documents/script_admb/",filename, sep=""))





